///@file SobolSequence.cpp
///@brief Implementation of the class representing a multidimensionnal Sobol Séquence
///@author Arnaud Duval
///@version 0.0
///@date 2010/10/27

#include "SobolSequence.h"
#include "ProbUtils.h"
#include <fstream>
#include <cstdlib>
#include <cmath>

///@brief Default constructor
SobolSequence::SobolSequence()
{
	D = 0;
	N = 0;
}

///@brief Constructor with given size and dimension
///@param d dimension
///@param n number of points
SobolSequence::SobolSequence(unsigned int d, unsigned int n)
{
	N = n;
	D = d;
	///@brief Reading file that contains the direction (obtained from Joe & Kuo)
	///@bug It would be better to work with relative path for this file
	std::ifstream infile("/home/arnaud/bergamote/data/joe-kuo-6.21000", std::ios::in);
	if(!infile)
	{
		///@warning A gestion of exception would be better
		std::cerr << "Unable to open file containing Sobol direction, exiting...\n";
		exit(1);
	}

	char buffer[1000];
	infile.getline(buffer,1000, '\n');

	///@brief Maximum number of necessary bits
	unsigned L = (unsigned)ceil(log((double)N)/log(2.));

	unsigned *C = new unsigned[N];
	C[0] = 1;
	for(unsigned i = 1 ; i < N ; i++)
	{
		C[i] = 1;
		unsigned value= i;
		while (value & 1)
		{
			value >>= 1;
			C[i]++;
		}
	}

	pts = new double * [N];
	for(unsigned i = 0 ; i < N ; i++)
		pts[i] = new double[D];
	for(unsigned j = 0 ; j < D ; j++)
		pts[0][j] = 0;

	///@brief Computation of first direction
	unsigned *V = new unsigned[L+1];
	for(unsigned i = 1 ; i <= L ; i++)
		V[i] = 1 << (32-i);

	unsigned *X = new unsigned[N];
	X[0] = 0;
	for(unsigned i = 1 ; i < N ; i++)
	{
		X[i] = X[i-1] ^ V[C[i-1]];
		pts[i][0] = (double) X[i]/pow(2.0,32);
	}

	///@brief Cleaning
	delete[] V;
	delete[] X;

	///@brief Compute other dimensions
	for(unsigned j=1 ; j < D ; j++)
	{
		///@brief reading parameters from file
		unsigned d, s;
		unsigned a;
		infile >> d >> s >> a;
		unsigned *m = new unsigned[s+1];
		for(unsigned i=1 ; i<=s ; i++)
			infile >> m[i];

		unsigned *V = new unsigned[L+1];
		if(L <= s)
		{
			for(unsigned i=1 ; i <= L ; i++)
				V[i] = m[i] << (32-i);
		}
		else
		{
			for(unsigned i=1 ; i <= s ; i++)
				V[i] = m[i] << (32 - i);
			for(unsigned i=s+1 ; i <=L ; i++)
			{
				V[i] = V[i-s] ^ (V[i-s] >> s);
				for(unsigned k=1; k<s ; k++)
					V[i] ^= (((a >> (s-1-k)) & 1) * V[i-k]);
			}
		}

		unsigned* X = new unsigned[N];
		X[0] = 0;
		for(unsigned i = 1 ; i < N ; i++)
		{
			X[i] = X[i-1] ^ V[C[i-1]];
			pts[i][j] = (double)X[i]/pow(2., 32);
		}
		///@brief Cleaning
		delete[] m;
		delete[] V;
		delete[] X;
	}
	delete[] C;
}


///@brief Destructor
SobolSequence::~SobolSequence()
{
	delete[] pts;
}

///@brief overload of [] operator : direct access to a realisation
///@param i index of realisation to return
double* SobolSequence::operator[] (unsigned int i)
{
	assert(i > 0);
	assert(i <= N);
	return pts[i-1];
}

///@brief Transform a variable to Gauss distribution
///@param j index of variable to transform
///@param av average value for Gauss distribution
///@param dev Standard deviation for Gauss distribution
void SobolSequence::TransformGauss(unsigned int j, double av, double dev)
{
	assert(j > 0);
	assert(j <= D);
	for(unsigned int i = 0 ; i < N ; i++)
	{
		pts[i][j-1] = RecFgauss(Fu(pts[i][j-1],0.,1.), av, dev);
	}
}

///@brief Transform a variable to a uniform distribution
///@param j index of variable to transform
///@param a lower bound for uniform distribution
///@param b upper bound for uniform distribution
void SobolSequence::TransformUniform(unsigned int j, double a, double b)
{
	assert(a <= b);
	assert(j > 0);
	assert(j <= D);
	for(unsigned int i = 0 ; i < N ; i++)
	{
		pts[i][j-1] = Recu(Fu(pts[i][j-1],0.,1.),a,b);
	}
}

///@brief Return number of points
int SobolSequence::Nbpts()
{
	return N;
}

///@brief Return dimension
int SobolSequence::Dim()
{
	return D;
}

///@brief Overload of << operator
///@param s out stream
///@param a Sobol Sequence to put inthe out stream
std::ostream& operator << (std::ostream& s, const SobolSequence& a)
{
	for(unsigned int i = 0 ; i < a.N ; i++)
	{
		for(unsigned int j = 0 ; j < a.D ; j++)
			s << a.pts[i][j] << "\t";
		s << "\n";
	}
	return s;
}


